Exact diagonalization of the one dimensional Bose-Hubbard model 

with local 3-body interactions 



(N 

o 

(N 

on: 



o3 



i 

o 
o 



> 

On 1 

(N 
O 
(N 



X 



Tomasz Sowinski 

Institute of Physics of the Polish Academy of Sciences, Al. Lotnikow 32/46, 02-668 Warsaw, Poland 

(Dated: February 10, 2012) 

In this article the extended Bose-Hubbard model with local two- and three- body interactions is 
studied by exact diagonalization approach. Shapes of the first two insulating lobes are discussed and 
values of the critical tunneling for which the insulating phase loses stability for different three-body 
interactions are predicted. 
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Experimental progress on trapping and manipulating 
ultracold atoms confined in optical lattices has opened 
new perspectives for controlling many-body states of dif- 
ferent quantum systems 0, Q • In the simplest case, when 
atoms interact mainly via two-body contact interactions, 
such systems are described in the context of the Bose- 
Hubbard (BH) model characterized by one free parame- 
ter only - the ratio of tunneling amplitude J to the on- 
site two-body contact interaction strength U. For such 
systems the quantum phase transition from superfluid 
(SF) to Mott insulator (MI) has been predicted and ob- 
served • Depending on physical realization there are 
many different extensions of standard BH models pre- 
dicting existence of novel quantum phases [5-9]. One 
of possible extensions of the BH originates in taking into 
account not only two- but also higher many-body interac- 
tions. It is highly nonintuitive how many-body terms can 
appear in a description of systems in which all mutual in- 
teractions are purely two-body. One should note however 
that direct interactions between particles are introduced 
as an effective description in the limit of low energy. 
Therefore it is theoretically possible that in some partic- 
ular experimental scenarios higher many-body terms can 
significantly change properties of the system. One of such 
scenarios was proposed for a polar molecules in optical 
lattices. As it is shown in [wi pLlj l it is possible to con- 
trol two- and three- body interactions between molecules 
independently. The other source of many-body terms in 
the presence of the optical lattice potential comes from 
the fact that mutual interactions can excite atoms to the 
higher orbital states and modify the shape of the single- 
particle wave functions. This phenomenon can be effec- 
tively taken into account in the model by adding three- 
body interaction terms. In the case of repulsive forces 
a single atom wave-functions are essentially broadened 
and therefore effective three-body terms become attrac- 
tive [H, 03 . 

In this article the ground state of the one dimensional 
extended Bose-Hubbard model describing atoms inter- 
acting locally via two- and three- body interactions is 
studied. Such a model was previously discussed in the 
mean field regime QJ, [H|. Additionally the first lobe 
of the insulating phase, i.e. when an average number of 
particles per lattice site p = 1, was recently studied with 



density matrix renormalization group (DMRG) methods 
fl6| . Here we investigate physically more interesting re- 
gion of the phase diagram when three-body interactions 
direct influence on the properties of the ground state, i.e. 
when p = 2. It is clear that in such a case, in contrast 
to the p = 1 case, any tunneling process which leads 
to destruction of the insulating phase has to concur not 
only with two-body interactions but also with three-body 
ones. Moreover, in all previous analysis only repulsive 
three-body terms were taken into account. In the light 
of recent experiments [lH it is clear that attractive case 
is probably more realistic and therefore I study also this 
case here. 

From now on I assume that higher orbital interactions 
can be effectively taken into account by adding a three- 
body interaction term. Therefore I restrict myself to the 
lowest band of the optical lattice and then the system of 
N bosons confined in a chain of L lattice sites with pe- 
riodic boundary conditions is described by the following 
extended Bose-Hubbard Hamiltonian 

l u L 

H = - J ^(ajfii+i + h.c.) + — y^nj(nj - 1) 
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The bosonic operator fij (aj) annihilates (creates) parti- 
cle at site i, J is the nearest neighbor hopping amplitude, 
and U is the on-site two-body repulsion energy between 
bosons. Additional term proportional to W describes the 
on-site three-body interaction. For convenience I intro- 
duce particle number operators hi — ajfij and average 
occupation parameter p = N / L where N is a total num- 
ber of particles in the system. Note, that the additional 
three-body term can be viewed as a first occupation- 
dependent correction to the on-site two-body interaction 
in the standard Bose-Hubbard Hamiltonian 



W 

U(n) = U+— (n-2) + 
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Occupation-dependence comes from the fact that for 
higher filling of the lattice site the lowest Wannier func- 
tion does not properly describe the localized orbital. Sim- 
ilarly, the ground state of a harmonic oscillator wave 
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function does not resemble the ground state solution of 
the Gross-Pitaevskii (GP) equation for large condensate 
fraction in the harmonic trap. Then instead of expand- 
ing the ground-state solution in the basis of trap poten- 
tial one can find the effective occupation-dependent fre- 
quency for which the harmonic oscillator ground-state 
mimics the GP solution. In the case of BH systems 
one can still use the lowest Wannier state but with 
occupation-dependent interactions (l2l. Ilijj. 

The Hamiltonian (TTJ) is studied for repulsive three- 
body interactions (W > 0) as well as for attractive one 
(W < 0). In the second case, it is necessary to take into 
account also a weak four-body repulsive interactions to 
prevent system against collapsing. 

In this brief report the diagram of quantum phases 
in the presence of the three-body interactions is stud- 
ied theoretically. The numerical approach is based on 
the exact diagonalization of the Hamiltonian (TTJ) and fol- 
lows the approach discussed first time by Elesin et al. in 
pjj adopted later to more sophisticated numerical tech- 
niques. The method draws on the observation that for 
infinite size of the system MI phase, in contrast to the 
SF, has non zero energy gap A for adding/subtracting 
particle to/from the system. Therefore, in principle one 
can determine insulating lobes in the phase diagram by 
finding ground state energies of systems differing by one 
particle. It was shown by Fisher et al. [18J that for 
models considering on-site interactions only the insulat- 
ing phase can occur only for integer p. In such a case 
the MI phase can be stable for nonzero hopping ampli- 
tude J. It can be easily shown with direct arguments 
that in the limit J — > the energy gap of the first lobe 
(p = 1) is not affected by three-body terms and it is 
equal to A = U. In contrast, the second lobe (p = 2) 
is highly influenced by three-body interactions and its 
width is equal to A = U + W for J = 0. To find en- 
ergy gap A for infinite system for given parameters U, 
W , and nonzero tunneling amplitude J we follow a sim- 
ple procedure. First, we perform exact diagonalization 
of the Hamiltonian for finite size L and fixed number of 
particles N. In this way we find the ground state of the 
system |G), as well as the ground state energy E(L,N). 
Then we define upper (p+) and lower limit of the 

insulating phase as follows 



p+{p,L) = E(L, P L^ 
l x_{p,L)=E{L,pL) 



l)-E(L,pL), 
-E(L,pL-l) 



(3a) 
(3b) 



These quantities highly depend on the size of the lattice. 
Since we are interested rather in their values in the limit 
of infinite system, therefore we diagonalize Hamiltonian 
for different lattice sizes (in our calculations we take L = 
5, ... 8). Then we plot p + and p_ as a functions of 
and we extrapolate the data to the point 1/L = 0. In 
this limit we define energy gap of insulating phase A as a 
difference A = p + — /i_ . To show how this method works 
in practice we present in Fig. Q] an example for W — U 
and J/U = 0.2 in two cases p = 1 and p = 2. As it is seen 
both mu's depend almost linearly on 1/L and therefore 
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7.16 x 10" 


2 2.85 x 10" 1 
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7.10 x 10" 


2 2.75 x 10" 1 



TABLE I: Energy gap of insulating phase in the limit of in- 
finite system size obtained for different fitting procedures. 
Presented values are calculated for specific parameters of the 
Hamiltonian (V = U and J/U — 0.2) discussed as an exam- 
ple. 




FIG. 1: The upper p+ and lower /i_ limit of the insulating 
phase as a function of the inverse of the system size 1 / L for 
two densities p — 1 (left panel) and p = 2 right panel (Hamil- 
tonian parameters are: V = U and J/U = 0.2). The solid 
lines is a linear fit to the numerical data points. Extrapola- 
tion to infinite system size 1/L — > gives energy gap A of the 
insulating phase in the limit of infinite system. 



linear data extrapolation gives quite good predictions. 
To test an accuracy of obtained energy gap A we fit also 
higher order polynomials in 1/L and compare resulting 
values. Predicted values of A for linear, quadratic, and 
cubic polynomial are presented in Table HI As it is seen a 
differences of obtained values of energy gap are equal to 
about 2% for p = 1 and about 10% for p = 2. Therefore, 
for further calculations a validity of linear dependence of 
limiting /i's and A on 1/L is assumed. 

With the procedure described above one can find insu- 
lating gap A for different parameters of the Hamiltonian 
and different densities (left panel in Figj2]). The transi- 
tion point from MI to the SF phase occurs when A be- 
comes equal to zero. Technically it is very hard to adopt 
this definition directly since A has some numerical uncer- 
tainty. In this paper it is assumed that this uncertainty 
comes mainly from the extrapolating procedure and it is 
equal to 5 x 10~ 3 C/. In a consequence we estimate a po- 
sition of the transition point as a hopping amplitude for 
which insulating gap becomes smaller than this uncer- 
tainty. These points are marked with open circles in Fig. 

The phase diagrams in the right panel were obtained 
by plotting p± as a functions of tunneling amplitude J. 
Additionally, in the background the density plot of the 
correlation function $ = (G\a\ai+i \G)/N is shown. The 
meaning of this quantity is given below. 

For W = we recover well known phase diagram 
for standard Bose-Hubbard model (upper panel in Fig. 
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[5]). Transition points for the first and the second MI 
lobes occur for J c ss 0.29J7 and J c « 0.17C/ respectively. 
These values are in good agreement with previous cal- 
culations [l9l [20l| . For positive (negative) W we observe 
that transition point for first MI lobe moves slightly to- 
wards larger (smaller) hopping amplitudes J and smaller 
(larger) chemical potential fi. Additionally, energy gap 
A in the limit of vanishing tunneling is equal to U for 
any W . It is quite obvious since in this limit every lat- 
tice site is occupied exactly by one boson and therefore 
three-body interactions can be completely neglected. In 
this way we confirm observations based on DMRG meth- 
ods |lq |. The situation is quite different for second in- 
sulating lobe. In this case, for repulsive three-body in- 
teractions (W > 0), the SF phase occurs for larger tun- 
neling amplitudes. It comes from the fact that tunnel- 
ing which destroys insulating phase has to overcome the 
energy barrier increased by W. In opposite, for attrac- 
tive three-body interactions (W < 0) destroying of the 
MI phase should be much simpler since then a three- 
body configuration is less bounded. It is very simple to 
show that in the limit of vanishing J the energy gap A 
(width of the MI lobe in the phase diagram in the fi di- 
rection) should be equal to U + W. It means that in the 
limit W — » — U a vanishing of the second insulating lobe 
should be observed. Numerical results fully confirm all 
these predictions (Fig. [5]). In Fig. [3] we show the depen- 
dence of critical hopping amplitude on the strength of 
the three-body interactions. As was suspected the crit- 
ical tunneling J c is much more sensitive to three-body 
interactions for second lobe than for the first one. 

Exact diagonalization method gives a quite nice con- 
firmation that MI phase can occur only for integer fill- 
ings. One can calculate the correlation function $ = 
(G\a i+1 ai\G) / N which describes hopping of particles be- 
tween neighboring sites and it might be viewed as number 
conserving analogous of the mean field. Due to the trans- 
lational invariance this quantity does not depend on the 
index i. From numerical calculations for finite lattices 
we find that for vanishing tunneling amplitudes J this 
quantity drops exactly to zero only for integer fillings. In 
other cases it remains nonzero even in the limit J — > 0. 
Moreover, the critical tunneling J c determined previously 
from the energy gap arguments occurs when $ (for inte- 
ger filling) starts to decay much faster than for surround- 
ing non integer ones. In Fig. @] values of $ for different 
fillings and three-body interactions are presented (num- 
ber of lattice sites L = 8). Additionally, in Fig. [5] the 
density plot of $ is visualized in the background of ap- 
propriate phase diagrams. Darker areas correspond to 
smaller values of $. In these cases 5> is calculated by ex- 
act diagonalization of the grand canonical Hamiltonian 
H — fiN in truncated basis with L = 8. The insulating 
lobes predicted this way match the lobes determined by 
energy gap boundaries. 

Finally, let me note that standard one dimensional 
Bose-Hubbard model has totally different properties than 
predicted by mean field approximation [21|. This fact 




FIG. 2: Properties of the system for different values of the 
three-body interaction strengths: W — (top), W/U = —0.5 
(middle) and W/U = 1 (bottom). On the left I present the 
energy gap A of the first (p = 1) and the second (p = 2) MI 
lobe as a function of hopping amplitude J. On the right the 
phase diagrams of considered systems are presented. Dotted 
lines represent values of p- and fi+ for infinite system ob- 
tained by extrapolation from finite system properties. Open 
circles are markers of the transition points defined as a points 
where A is equal to zero at the edge of significance. In the 
background of each phase diagram the density plots of the 
correlation function $ are visualized. As was expected for 
p = 1 insulating phase do not change significantly, but for 
p — 2 size of the MI phase crucially depends on the three- 
body interactions parameter. 
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FIG. 3: Values of the critical tunneling as a function of the 
three-body interaction parameter W for p — 1 (left panel) and 
p = 2 (right panel). In the first case, position of the critical 
point in the phase diagram remains almost unchanged as was 
previously predicted pq |. In the second case position of the 
critical point highly depends on the three-body interactions. 



originates in the observation that at the tip of the Mott 
lobe the one dimensional BH model belongs to the uni- 
versality class of the XY spin model. It means that the 
phase transition from MI to SF phase is of the Kosterlitz- 
Thouless (KT) type. It is characterized by the exponen- 
tial decaying of the correlation length [22]. In other words 
the energy gap A decays exponentially in the neighbor- 
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FIG. 4: The correlation function $ = (G|oJ +1 Oi|G)/JV as a 
function of tunneling amplitude J for different Hamiltonian 
parameters and fillings. It is seen that for integer fillings (bold 
solid lines), in contrast to other non integer ones (dashed and 
thick solid lines), function $ decays to zero in the limit J — > 0. 
The transition tunneling J c predicted by other method occurs 
when integer filling <£> detaches from the others. This behavior 
is universal for repulsive (W > 0) as well as for attractive 
(W < 0) three-body interactions. Note that for plots legibility 
a nonlinear scaling is used. 



hood of the critical point, i.e. ln(A/J7) is a linear function 
of l/\/l — J I J c for J < J c . Numerical results based on 
exact diagonalization show that these behavior does not 
depend on the three-body interactions and KT transition 
occurs in the system also in the case when three-body in- 
teractions (attractive as well as repulsive) are taken into 
account. 

To conclude, in this brief report the extended Bose- 
Hubbard model with additional local three-body inter- 
actions was studied. To find a critical behavior of the 
system the exact diagonalization of the Hamiltonian was 
used. In this way, the phase diagrams for repulsive and 
attractive three-body interactions were obtained. It was 
shown that shape of the second insulating lobe (for two 
bosons per lattice site in average), in contrast to the first 
one, crucially depends on the three body interactions. 
Additionally, it was checked that local three-body inter- 
actions do not change the characteristic properties of the 
phase transition. 
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